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COMPUTER ANALYSIS OF MULTICIRCUIT 


SHELLS OF REVOLUTION BY THE FIELD METHOD 
By Gerald A. Cohen 

Structures Research Associates, Laguna Beach, CA 


SUMMARY 


The "field method", presented previously for the solution of even- 
order linear boundary-value problems defined on one-dimensional open 
branch domains (i.e., trees), is extended to boundary-value problems 
defined on one-dimensional domains containing circuits (i.e., multicircuit 
graphs) . This method converts the boundary-value problem into two 
successive numerically stable initial-value problems, which may be solved 
by standard forward integration techniques. In addition, a new method for 
the treatment of singular boundary conditions (i.e., kinematic constraints) 
is presented. This method, which amounts to a (partial) interchange of 
the roles of "force" and "displacement" variables, is problem independent 
with respect to both accuracy and speed of execution. 

The method has been implemented in a computer program to calculate 
the static response of ring-stiffened orthotropic multicircuit shells of 
revolution to asymmetric loads. Solutions are presented for sample 
problems which illustrate the accuracy and efficiency of the method. 


INTRODUCTION 


This paper is a sequel to an earlier paper (ref. 1) in which it is 
shown how linear boundary-value problems defined on one-dimensional 
branched domains without circuits (i.e., trees) may be formulated and 
solved as two successive initial-value problems. In contrast to the more 
conventional superposition methods for such problems, the integration of 
the resulting initial-value problems is numerically stable, thereby 
eliminating the need for artificial subdivision of the tree into a number 
of "suitably small" subintervals. This method of solution has been termed 
the "field method"; however, the author has since become aware of a 
rapidly expanding body of related literature under the general heading of 
"invariant imbedding" (ref. 2). There are several variations of the 
invariant imbedding method for solving boundary-value problems, and the 
field method (variously known as the method of sweeps, the method of 
factorization, or the method of generalized Riccati transformations) may 
be considered to be one of them (ref. 3). Although the alternate forms of 
invariant imbedding eliminate the second (backward) integration of the 



field method, they suffer from the need of having to choose a priori the 
set of response output (or storage) points. In the field method , since 
the response functions are obtained directly from the backward integration, 
their storage points may be chosen automatically so as to obtain a 
uniformly dense description of these functions. Furthermore, in many 
cases the extra work required in the forward integration of the alternate 
methods more than compensates for the elimination of the backward 
integration. 

In reference 1 the field method is applied to the solution of the 
equations of static response of open branch ring-stiffened shells of 
revolution. Because the field method not only eliminates the "long 
subinterval" (i.e., numerical instability) problem of superposition 
methods, but is shown, by example, to execute significantly faster, it 
appears to be a desirable method for shell problems. The purpose of the 
present paper is twofold: 1) to generalize the field method to boundary- 
value problems defined on one-dimensional domains which contain circuits, 
and 2) to improve the treatment of singular boundary conditions (i.e., 
kinematic constraints) given in reference 1, in which the numerical 
solution depends on a scaling matrix, an effective choice of which is 
problem dependent. 


SYMBOLS 


a,b,c,d matrix coefficients of linear differential equations [eqs. (1)] 

B,D matrix coefficients of linear boundary conditions [eq. (2)] 

D,L effective values of D and L at interior vertices [eqs. (39) 

and (41) ] 

D,L effective values of D and L at initial vertex of a 

continuation arc [eqs. (53) and (59)] 

f,g inhomogeneous vectors of linear differential equations 

[eqs. (1)] 

G.P. generic point of a closed branch arc 

I identity matrix 

I.P. "initial point" of a closed branch arc (see p. 6) 

K stiffness matrix for part of disconnected graph 

between the I.P. and G.P. 

L inhomogeneous vector of linear boundary conditions 

[eq. (2)] 


Mi 


meridional stress couple 
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n circumferential harmonic number 

P,Q,S shell forces per unit of circumferential length in 

axial, radial, and circumferential directions 

p half -order of boundary-value problem 

R -2 circumferential radius of curvature 

r small circle radius 

s arc distance 

s ,<p,z meridional, circumferential, and normal coordinates 

of shell reference surface 

t wall thickness 

u, w field function matrices [eqs. (6a) and (7)] 

v, <S,y»A additional field function matrices for closed branch 

arcs [eqs. (6)] 

x, y axial and radial coordinates 

y, z generalized force and displacement response vectors 

a, 3 diagonal matrices with elements equal to zero or one 

[eqs. (20) and (21)] 

A net change across a vertex [eq. (3)]; also effective 

value of D at I.P. of arcs in a split circuit [eq. (51d) ] 

k B^D 

y,ViQ,IJ 2 inhomogeneous vector and matrix coefficients in equation 

for z at I.P. of arcs in a split circuit [eqs. (51)] 

£,ri,v shell displacements in axial, radial, and circumferential 

directions 

u effective value of v at initial vertex of a continuation 

arc [eq. (53b)] 

X meridional rotation 

Subscripts : 

A,E value at vertex P on arcs A and E of figure 6(b) 

C value at vertex P on circuits C of figure 6(b) 
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P value at vertex P of figure 6(b) 


0 

value at I.P. of a 

closed branch arc 

1,2 

value at vertices 

Pj and P 2 of figure 6(a) 

Superscripts : 


T 

matrix transpose 


( ) T 

d( )/ds 


( ) + 

value at vertex on 

exiting arc 

( )" 

value at vertex on 

entering arc 

(“) 

modified variable 

for singular arcs 


ANALYTICAL FORMULATION 


As in reference 1, in order to formulate the field method in a 
general context, it will be convenient to employ some elementary concepts 
from the theory of geometric graphs. Figure 1 shows the reference 
meridian of a hypothetical shell of revolution. The heavy dotted points 
depict boundaries of the shell, which are defined as one of the following 
types of points: 

(1) branch points 

(2) branch extremities 

(3) ring or ring load points 

(4) shell property or load discontinuity points 

The boundaries thus divide the meridian into a number of sub intervals . 

In contrast to other numerical integration methods, no additional 
(artificial) boundaries are required simply to reduce subinterval length. 

A geometric graph (ref. 4) is defined as a set of points, called 
vertices, and a set of non-self-intersecting curves, called arcs, 
satisfying the following requirements: 

(1) Each closed arc contains precisely one vertex. 

(2) Each open arc contains precisely two vertices, viz. its 
end points. 

(3) The arcs have no common points, except for the vertices. 

It is clear from figure 1 that if we identify the boundaries as vertices 
and the subintervals as arcs, the reference meridian of a shell of 
revolution is nothing more than a geometric graph. 
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The following graph terminology will be used: 

(1) Chain - a continuous sequence of arcs from an initial vertex 
to a terminal vertex. A non-self-intersecting chain is said 
to be simple. 

(2) Circuit - a simple chain whose initial and terminal vertices 
coincide. 

(3) Connected graph - a graph for which every pair of vertices is 
joined by at least one chain. 

(4) Tree - a connected graph which contains no circuits. 

Whereas the analysis of reference 1 is limited to boundary-value 
problems defined on trees, the present analysis treats a much broader 
class of connected graphs which may contain multiple circuits. In fact, 
the sole limitation on the generality of graphs treated is contained in 
the following statement. If two vertices of the graph are joined by 
three or more distinct simple chains, then all circuits containing one 
of these vertices also contain the other. The reason for this limitation 
as well as its implications, are discussed later in the paper (see p. 16) 
For each arc of the graph, the arc distance s will be used as the 
independent variable. Let us assume that the arcs have been ordered and 
oriented (with respect to the direction of increasing s) , but defer for 
the time being the manner in which this was done. 


Definition of Boundary-Value Problem 

A system of ordinary differential equations of order 2p may always 
be written as a system of 2p first-order equations. If we group one-half 
of the dependent variables in the p x 1 matrix y and the other half in 
the p x X matrix z, the system of first-order equations may be written, 
for linear systems, as two matrix differential equations, viz. 

y f + ay + bz = f (la) 

z’ + cy + dz = g (lb) 

where prime denotes differentiation with respect to s; a, b, c, and d are 
p x p matrix functions of s; and f and g are p x 1 matrix functions of s. 

Equations (1) are defined at every interior point of each arc of the 
graph. They are supplemented by linear boundary conditions, defined at 
the vertices of the graph, of the form 

BAy + Dz = L (2) 

where B and D are p x p matrices, L is a p x 1 matrix, and 

Ay = I y + - I y" (3) 

Here, y and y represent the values of y at the vertex on exiting (s 
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increasing away from the vertex) and entering (s increasing towards the 
vertex) arcs, respectively. As implied by the form of eq. (2), it is 
assumed that z is continuous at vertices.* A vertex and its boundary 
condition are said to be singular if the matrix B is singular; otherwise 
they are called regular. t 


The majority of boundary-value problems in mechanics are self-adjoint. 
This property is synonymous with being derivable from a variational 
principle. Physically, this corresponds to systems which do not involve 
energy losses. For the system defined by eqs. (1) and (2), the conditions 
of self-adjointness may be shown to be 


b 



K = 


T 


< 



(4a) 

(4b) 


where 


k = B“ ] D 


(5) 


and the superscript T denotes matrix transpose. Although condition (4b) 
strictly applies only to regular vertices, it may be used also for 
singular vertices if a singular vertex is viewed as a limiting case of 
a sequence of regular vertices. Thus, for self-adjoint problems a 
singular vertex is the limit of a sequence of regular vertices, for 
each of which eq. (4b) holds true. 


Field Relations 

An open branch arc is characterized by the fact that a cut at any of 
its points disconnects the graph into two separate parts. Conversely, a 
(single) cut of a closed branch arc (not to be confused with a closed arc , 
defined on p. 4), which, by definition, lies in a circuit, does not 
disconnect the graph. However, in order to develop field relations , it 
is necessary to visualize disconnecting the graph by a cut at a generic 
point. For this purpose, the concept of an Initial point , denoted hence- 
forth by I.P., of a closed branch arc is introduced. For each closed 
branch arc, the I.P. is chosen so that cuts at the I.P. and a generic 
point, denoted henceforth by G.P., of the arc disconnect the graph into 
two separate parts. Note that the I.P. of a given closed branch arc is 
located in a circuit containing the given arc but not necessarily in that 
arc itself, and in fact several closed branch arcs lying in the same 
circuit may have the same I.P. (fig. 2). 


*It is assumed on physical grounds that one-half of the dependent variables 
are continuous at vertices, and these shall be grouped in the vector z. 

The vector y may be viewed as a generalized force vector corresponding to 
the generalized displacement vector z, and -^Ay represents the net external 
force entering the Vertex. 

tA singular boundary condition implies a relationship between the components 
of the displacement vector z, i.e., kinematic constraint. 
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Since cuts at the I.P. and G.P. of a closed branch arc disconnect 
the graph Into two separate parts, it is clear that the values z - zq 
at the I.P. and z - z at the G.P. can serve as boundary conditions which, 
along with the differential equations (1) and the boundary conditions (2) 
over one of the parts , determine y and z over that part independently of 
the corresponding data over the remaining part.* In particular, the 
values zq and z determine the values of y at the cuts, which for linear 
problems is expressed by the matrix relations 


y = uz + vz o + w 

(6a) 

4 1 

-y 0 = 6z + yzo + X 

(6b) 


where u, v, 6, and y are p x p matrices and w and A are p x 1 matrices. 
Equations (6) are the field relations for closed branch arcs, which 
generalize eq. (6) of reference 1 for open branch arcs, viz. 

y = uz + w (7) 

■j- 

In eq. (6b), -y q represents the force vector at the I.P. on an exiting 
closed branch arc, which shall be called an initial closed branch copc, 

[All nonsub scripted variables in eqs. (6) are understood to be evaluated 
at the G.P.] Thus, by definition, the I.P. of an initial closed branch 
arc coincides with its initial vertex. In order to obtain an initial- 
value problem for the field functions u, v, w, 6, y, and A, closed branch 
arcs are ordered and oriented so that for every cut pair (at the I.P. and 
the G.P.) all of one part of the disconnected graph is completely 
described by s-values smaller, than that at the G.P. (fig. 2).+ 

Several properties of the field functions may be deduced on physical 
grounds. Because of the generally decaying nature of shell response to a 
local perturbation, it follows that as the G.P. recedes from the I.P., the 
coupling between response variables at these two points becomes weak. 

From eqs. (6), this is expressed by the statement that the matrices v and 
6 decay in magnitude, or symbolically 


v -»■ 0 

(8a) 

S + 0 

(8b) 


as the distance between the G.P. and the I.P. becomes large. From eq. 


*It is assumed that giving z q and z at the cut pair determines a unique 
solution, i.e., there exists no nontrivial solution of the homogeneous 
forms of eqs. (1) and (2) satisfying the cut conditions zq = z * 0, as 
is the case for problems derivable from the minimization of a positive- 
definite functional. Otherwise, a transformation of variables, such as 
used on singular arcs (see p. 11), may be necessary. 

tin so doing, it will be assumed that the points of the fully described 
part are jput into one-to-one correspondence with s-values between the 
values of> s at the I.P. and at the G.P. 
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(8b) it follows that eq. (6b) cannot be used for a direct calculation of 
z, once yj and zo (and the field functions) are known. For solution of 
eq. (6b) for z gives 

z = “6 _1 (yJ + yz 0 + A) (9) 

Since z does not grow large with 6 -1 (and since y and A do not decay and 
y"o and zq are constants), it follows that at points remote from the I.P. 
the calculation of z from eq. (9) would involve the small difference of 
large numbers.* 


For self-adjoint problems, to which our attention is henceforth 
confined, the matrix K defined by 



( 10 ) 


represents a stiffness matrix for the part of the disconnected graph 
which is fully described by s-values smaller than s at the G.P. Since 
such a stiffness matrix is known to be symmetric, its submatrices 
satisfy the following reciprocal relations 

T T , T /11N 

u = u , Y = Y, S=v (11) 


Thus, for example, for eighth-order boundary-value problems, there are 
44 independent (scalar) field functions contained in eqs. (6). However, 
as shown below, only 30 of them (u, v, and w) need be stored. 


Plan of Field Method 

Before proceeding to derive the differential equations and initial 
conditions for the field functions, in order to put the discussion in 
context, a preview of the field method for a graph with circuits is 
presented. The calculation procedure consists of the following three 
basic steps (fig. 3). 

I A forward integration is performed over the graph for u and w on 


*Not only do each of the elements of v and 6 approach zero but, in general, 
the matrices themselves become numerically singular. This follows from 
the fact that the columns of v and the rows of 5 are solutions of a linear 
homogeneous system of differential equations [see eqs. (14c) and (15)]. 

Each of the columns of v and the rows of 6 is therefore a linear combina- 
tion of a fundamental set of (decaying) solution vectors. As the G.P. 
recedes from the I.P., the fundamental solution vector which decays the 
slowest will dominate the other components so that, regardless of their 
initial values, the columns of v and the rows of 6 will lose their 
numerical independence. Thus, the inversion of 6 will, by itself, incur 
increasing round-off error, and hence any scheme of direct calculation of 
z using 6 _1 will fail. 
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II 


open branch arcs and u, w, v, y, and A on closed branch arcs. On open 
branch arcs u and w are stored, and on closed branch arcs u, w, and v are 
stored. 

II Upon reaching the final vertex, in principle, two cases must be 
distinguished. If the final arc is an open branch arc or a closed branch 
arc whose I.P. coincides with its terminal vertex, the field relations and 
boundary conditions at the final vertex are solved simultaneously to give 
the value of z, say Z]_, at the final vertex. If, on the other hand, the 
final arc is a closed branch arc whose I.P. is distinct from its terminal 
vertex, then the field relations and boundary conditions at the final 
vertex and the I.P. are solved simultaneously to give z^ and the value of 
z at the I.P. 

III Using zi as a final value, and z-continuity at interior vertices, 
a backward integration of the equation [cf. eqs. (lb) and (7)] 

z* + (cu + d)z = g - cw (12) 

is performed over the graph. Before integrating eq. (12) on a closed 
branch arc, w on that arc is replaced by w + vzq where zq is the value 
of z at its I.P. [cf. eqs. (6a) and (7)]. As shown later, zq for each 
closed branch arc can be calculated prior to the backward integration on 
that arc. When a response storage point is reached, y is calculated from 
the field relation (7). 


Differential Equations for Field Functions 

The differential equations for u, v, and w are obtained by the same 
procedure used in reference 1 for u and w on open branch arcs. That is, 
eq. (6a) is differentiated with respect to s, the differential equations 
(1) are used to eliminate y' and z r , and eq. (6a) is used again to 
eliminate y. This gives 

(u’ - ucu +au-ud+b)z+ (v* - ucv + av)zg 

+ w T - ucw + aw + ug - f = 0 (13) 

Since eq. (13) must be satisfied identically with respect to z and zo 
(which depend on data at points of greater s, whereas u, v, and w do not), 
it follows that 

u f - ucu + au - ud + b = 0 (14a) 

w T - ucw + aw + ug - f = 0 (14b) 

v 1 - ucv + av — 0 (14c) 

Note that eqs. (14a, b) apply to both open and closed branch arcs, and that 
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eq. (14c) is simply the homogeneous form of eq. (14b).* 

In a similar manner, to obtain the differential equations for y and 
X, differentiate eq. (6b) with respect to s, and use eq. (lb) to 
eliminate z’ and eq. (6a) to eliminate y. This gives 

[<5 f - 5(cu + d)]z + (y ? - 5cv)zq + X' + 6 (g - cw) = 0 (15) 

Since eq. (15) is also an identity in z and zo, it follows from eqs. (15) 
and (11) that [the coefficient of z in eq. (15) being identically zero 
since by eqs. (4a) and (11) it is the transpose of the left side of eqs. 
(14c)] 

y f - v^cv = 0 (16a) 

X’ + v^(g - cw) = 0 (16b) 

These are remarkably simple differential equations, which may be integrated 
by quadrature after the integration for v and w is made. 


Initial Values of Field Functions 

In general, there are three different types of arcs for which initial 
values of the field functions must be obtained. These shall be called: 

(1) ordinary arcs, (2) initial closed branch arcs, and (3) continuation 
arcs . 


Ordinary arcs .- An ordinary arc is either: (1) an open branch arc 
whose initial vertex is incident with no closed branch arcs, or (2) a 
closed branch arc whose initial vertex is incident with precisely one 
other closed branch arc which precedes the ordinary arc and enters (i.e., 
is oriented towards) its initial vertex (fig. 4). For example, in figure 
2 arc 3 is an ordinary arc. Note that in the case of an ordinary closed 
branch arc, its I.P. is also the I.P. for the preceding closed branch arc. 

Regular arcs: An ordinary arc is called regular if its initial 
vertex is regular. The initial values of the field functions u, v, and w 
on a regular ordinary arc are obtained by substituting the field relation 
(6a) into the boundary condition (2) for its initial vertex to obtain 

(Au + B _1 D)z + Avz 0 + Aw - B _1 L = 0 (17) 

Since this is an identity in z and zo > and since there are no v- 
contributions from open branch arcs, eq. (17) yields the initial values 


*Henceforth, all results pertaining to u and w (or u and w) apply to both 
open and closed branch arcs. On open branch arcs, the closed branch 
field functions v, y, and X (or v, y, and X) may be taken to be zero and 
equations for them ignored. 
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u + = -B^D + l u~ 

(18a) 

w + = B“^L + £ w 

(18b) 

+ — 


V = V 

(18c) 


Moreover, evaluation of eq. (6b) on either side of the initial vertex, 
and subtraction of the results gives 


+ 

r = r 


x 


+ 


X 


(19a) 

(19b) 


Thus, all of the additional closed branch field functions are continuous 
at the initial vertex of an ordinary closed branch arc. 

Singular arcs: On an ordinary arc which is singular (|b| =0), eqs. 
(18a, b) break down. This problem is treated in reference 1 by forming 
on singular arcs a modified z-vector as a linear combination of z and y. 

In this relation, a scaling matrix, introduced as the coefficient of z, 
is required for dimensional homogeneity. Since this matrix must be 
determined empirically, its numerical effectiveness cannot be assured for 
all problems. In fact, the scaling matrix suggested in reference 1 
results in large round-off errors, as well as a slow forward integration, 
on singular arcs for shells of revolution with very high harmonic loading 
(viz. circumferential wave numbers n - 1000). For this reason, a universal 
transformation of variables for singular arcs was sought. 

In the case of an arc with initial clamping (i.e., the constraint 
z = 0) , the problem is resolved simply by an interchange of the roles of 
y and z. The generalization of this complete interchange of variables 
for an arbitrary singular arc, which may have at its initial vertex only 
partial kinematic constraint, is an interchange of some of the correspond- 
ing components of y and z. This is expressed by the transformation* 

y = ay + 3z (20a) 

z = ~3y + az (20b) 

where a and $ are diagonal p x p matrices with diagonal elements equal 
to zero or unity such that 

a + 6 = I (21) 


Equations (20) correspond to the interchange of the roles of the i-th 
components of the y and z vectors if the i-th diagonal element of $ is 
unity. As a rule, the i-th diagonal element of 3 is chosen to be unity 


*The minus sign in eq. (20b) is necessary in general to preserve the self- 
adjoint form of the differential equations (1) when written in terms of 
y and z [cf. eqs. (36) and (37)]. 
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if in order to make B nonsingular, it is necessary to replace the i-th 
column of B by the i-th column of D (see also p. 15). This is equivalent 
to interchanging components of y with corresponding components of z in 
sufficient number so that the coefficient matrix of y"*~ in the boundary 
equation (2) written in terms of y + and z + is nonsingular.* 



For a and 3 as defined, the following useful relations hold, 
aa = a 

(22a) 


33 = 3 

(22b) 


a3 = 3a = 0 

(22c) 

Using 

eqs. (21) and (22), eqs. (20) may be easily inverted to give 



y = ay - 3z 

(23a) 


z = 3y + az 

(23b) 

Thus, 

the inversion of eqs. (20) is effected simply by changing the 

sign 

of 3. 

Substitution of eqs. (23) into the field relations (6) gives 

the 

modified field relations 



y = uz + vzq + w 

(24a) 


-yt = 6z + yz 0 + A 

(24b) 

where 

u = (a - u3) _1 (ua +3) 

(25a) 


w = (a - u3) -1 w 

(25b) 


v = (a - u3) -1 v 

(25c) 

and 

6 = 5 (3u + a) 

(26a) 


y = y + 63v 

(26b) 


A = A + 63w 

(26c) 


Comparing the form of eqs. (20) and (24) with that of eqs (23) and (6), 
respectively, shows that eqs. (25) and (26) also are inverted simply by 
changing the sign of 6, i.e., barred and unbarred variables in eqs. (25) 
and (26) may be interchanged if 3 is replaced by -3. 


*That such an interchange is always possible is based on two premises: 

1) the boundary conditions (2) constitute four linearly independent 
equations, and 2) the boundary-value problem is self-adjoint. Note that 
the assumption of linear independence implies that the rank of the p x 2p 
matrix [B|D] is p, i.e., it contains a p x p nonsingular submatrix. 


12 



For self-adjoint problems, it is easily shown that the reciprocal 
relations (11) are satisfied by the modified field functions. To show 
that u is symmetric, start with the identity [ see eq. (21)] 


(a + 3)u = u(a + 3) (27) 

or au - u8 = ua - 8u (28) 

In view of eqs. (22), eq. (28) can be factored to 

(a - u8)(au + 3) = (ua + 3) (a - 8u) (29) 

or (au + 3) (a - 3u) -1 = (a - u3) _1 (ua + 3) (30) 


Because of the symmetry of u, a, and 3, the left-hand side of eq. (30) is 
the transpose of its right-hand side, which by eq. (25a) is u. Hence 

u T = u (31) 

To show that 6 is the transpose of v, use eqs. (11) and (31) to rewrite 


eq. (26a) in the form 

6 = v^(8u + a) = [ (a + u8)v]^ = v^ (32) 

where the last equality follows from the inverse of eq. (25c), viz. 

v = (a + u8) -1 v (33) 

As a consequence of eqs. (11), eqs. (26b, c) may be written as 

y = Y + v T 8v (34a) 

A = A + v T 8w (34b) 


and, as a consequence of eq. (32), barred and unbarred variables in eqs. 
(34) may be interchanged if 8 is replaced by -8. Thus from eq. (34a) and 
its inverse, one gets 

y - y = v 8v = v 8v - (v 8v) (35) 

which shows that the difference y ” Y is symmetric. Since, by eqs. (11), 

Y is symmetric, it follows that y is also symmetric. 

In terms of y and z, the differential equations (1) become 

y T + ay + bz = f (36a) 

z ' + cy + dz = g (36b) 


where 
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a = aaa + abB + Bca + BdB 

b = -aa3 + otba - 3c B + Bda 

c — -Baa - 3b 3 + aca + ad 3 

d = BaB - Bba - acB + ada 

f = af + 3g 
g = -Bf + ag 


( 37 ) 


Using eqs. (4a), it is easily shown from eqs. (37) that eqs. (4a) are 
satisfied also in terms of barred matrices, so that the transformation 
(20) preserves the self-adjoint property of the differential equations. 
Equations (24) and (36) can be used to derive differential equations for 
the modified field functions in exactly the same way as (6) and (1) were 
used to derive (14) and (16). Since eqs. (24) and (36) are identical in 
form to (6) and (1), respectively, the differential equations for the 
modified field functions are in the same form as eqs. (14) and (16) with 
all variables replaced by corresponding barred variables. 


Initial values of the modified field functions on singular ordinary 
arcs may be obtained from eqs. (18), (19), (25), and (34). As in reference 
1, we consider the singular vertex as the limit of a sequence of regular 


vertices, for each of which [cf. eq. (18a)] 

(a - u + 3) -1 = [a - (-B _1 D + £ u )B] -1 

= (Ba + D3) -1 B (38) 

where D = D - B £ u (39) 

Then, evaluation of eqs. (25) and (34) at the initial vertex of the 
ordinary arc gives, in view of eqs. (18) and (19), the initial values 

d + = -(Ba + D3) -1 (Da - BB) (40a) 

w + - (Ba + DB) - ^! (40b) 

v + = (Ba + D3) _1 Bv (40c) 

y + - -y + (v ^Bv"*" (40d) 

\ + = \~ + (v“) T Bw + (40e) 

where L = L + B £ w (41) 


Note that the transformation (20) includes as special cases the 
identity transformation 
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y = y, z = z (8 = 0) (42) 

used on regular arcs, and the complete interchange of variables 

y = 2 , z = -y (3=1) (43) 

used on arcs with initial clamping. Therefore, for the sake of unifying 
the method, all arcs could be integrated in terms of barred variables 
using the general equations (36), (37), and (40), depending on an 
appropriate choice for 3. However, because of the extra work involved in 
setting up the barred matrices of eqs. (37), it is advantageous in 
programming the method to retain the special transformations (42) and 
(43) as separate options.* Moreover, the complete interchange (43) is 
adequate to treat all singular arcs for which the specification of y + 
and, for closed branch arcs, zq uniquely determines z + . For problems 
derivable from the minimization of a positive-definite functional, this 
is the case, regardless of boundary conditions, for all closed branch 
arcs, and also for open branch arcs for which the initial vertex is an 
interior vertex if rigid body displacements are not admissible. For such 
problems, in practice, the only time the general (partial interchange) 
transformation (20) need be used is for singular open branch arcs whose 
initial vertex is a branch extremity or at whose initial vertex less than 
complete rigid body constraint is prescribed. t 

Initial closed branch arcs .- As discussed previously, an initial 
closed branch arc is characterized by the fact that its I.P. coincides 
with its initial vertex. Since the use of a minimum number of I. P.'s 
leads to the simplest analysis, in general an initial closed branch arc 
will be any closed branch arc for which no vertex preceding its initial 
vertex can serve as its I.P. A little reflection will show that a 
sufficient condition for a closed branch arc to be an initial closed 
branch arc is that its initial vertex is incident with another closed 
branch arc which lies in a common circuit and is either 1) an exiting 
arc or 2) an entering but not preceding arc. For example, in figure 2 
arcs 1, 2, 4, and 6 are initial closed branch arcs. 

Since the field relations (6) break down at the initial vertex of 


*As an example, for ortho tropic shells of revolution without property 
variation through the wall thickness, the calculation of the derivatives 
u' and w f [eqs. (14a, b)] consumes roughly 18 percent more time for the 
complete interchange (43) than for the identity transformation (42) and 
13 percent more time for the general transformation (20) than for the 
complete interchange, when the complete interchange and identity trans- 
formations are programmed as separate options. 

fActually, it suffices to have prescribed complete rigid body constraint 
at any vertex on the part of the disconnected graph fully described by 
s-values less than s at a cut on the singular arc. 
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an initial closed branch arc, such an arc is always treated as a singular 
arc, i.e., the integration on it is in terms of modified field functions. 
As for all singular closed branch arcs (see discussion above), the 
complete interchange of variables, eqs. (43), is used for initial closed 
branch arcs. The values of the field functions at the I.P. are obtained 
by inspection of the modified field relations (24) evaluated at the I.P. 
Substituting eqs. (43) into eqs. (24) and evaluating at the I.P. gives 


_ + 

zq m -u 0 yo + v 0 z 0 + w 0 
*yo = -^oyo + yo z o + x 0 


(44a) 

(44b) 


Since eqs. (44) must be satisfied identically with respect to yo and zq» 
it follows that 


u 0 = Yo = 0 

(45a) 

w 

II 

0 

(45b) 

w 0 - A 0 - 0 

(45c) 


Note that the additional initial value 6 0 = I implied by eq. (44b) is 
consistent with eq. (32) and is not used since, as noted earlier, for 
self-adjoint problems there is no need to calculate 6. 

Continuation arcs .- A continuation arc is defined simply as any arc 
not of the preceding two types. It is therefore either: 1) an open 
branch arc whose initial vertex is incident with a closed branch arc, or 
2 ) a closed branch arc whose initial vertex is incident with more than 
one other closed branch arc and which is not an initial closed branch arc. 
For example, in figure 2 arcs 5 and 7 are continuation arcs. 

In order to simplify the treatment of continuation arcs, a limitation 
of graph generality is introduced here. Henceforth our attention is 
restricted to graphs for which the following statement is true. If two 
vertices P\ and P 2 cere joined by three or more distinct simple chains , 
then all circuits containing P± (or P 2 ) also contain P 2 (or P\) . Examples 
of graphs admitted by this limitation, as well as more complex graphs not 
admitted, are shown in figure 5. (Note that the existence of additional 
open branches does not affect a graph's admissibility.) The treatment of 
such graphs reduces to the treatment of the elementary split circuit 
configuration shown in figure 6 (a). In this figure, arcs A and E may be 
open branch, closed branch, or even nonexistent arcs, and chains Cj and 
0-2 are representative of an arbitrary number of chains (i = l, 2 ,...,m) 
joining Pi and P 2 . According to this limitation, if A and (hence) E are 
closed branch arcs or C 3 exists, then no circuits, other than those 
depicted, containing Pi or P 2 exist. Additional open branch arcs incident 
with (and oriented towards) Pi or P 2 may exist and are accounted for in 
the analysis. 

The loop configuration shown in figure 6 (b) is the limiting case of 
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that in figure 6(a) when one of the chains approaches zero length, so 
that the vertices Pi and P 2 collapse into a single vertex P. It is 
necessary to include this limit configuration as a special case in the 
analysis since, as will be shown, the calculation of initial conditions 
for the continuation arc E in figure 6(a) entails ever-increasing round- 
off errors as Pj approaches P 2 . In figure 6(b) the circuit C is repre- 
sentative of an arbitrary number of such circuits containing the vertex P, 
and additional open branch arcs may enter P. For graphs consistent with 
the above-stated limitation, the sufficient condition for an initial 
closed branch arc given on page 15 is also a necessary condition, i.e., 
the initial vertex of an initial closed branch arc is incident with 
another closed branch arc which lies in a common circuit and is either 
1) an exiting arc [fig. 6(a)] or 2) an entering but not preceding arc 
[fig. 6(b)].* 

Split circuits: The analysis of the configuration of figure 6(a) has 
two objectives: 1) to determine the initial values of the field functions 
on the continuation arc E, and 2) to express z at the I.P. P^ in terms of 
z at P 2 in preparation for the backward integration of the chains C-^ (cf. 
step III, p. 9). The field relations available for this purpose are the 
following [cf. eqs. (6)]: 


Arc A evaluated at P^ 


yi = UiZi + VjZQ + W]_ 


(46a) 

+ - T - 

-yo = ( v i> z i + yi z 0 + 

X 1 

(46b) 

Chains evaluated at P 2 



y 2 = u 2 z 2 + v 2 zi + w 2 


(47a) 

+ - T - 

-yi = (v 2 ) z 2 + Y 2 z 1 + 

1 CM 
-< 

(47b) 

Arc E evaluated at P 2 



+ + , + , + 

y 2 = u 2 z 2 + v 2 Zq + w 2 


(48a) 

+ , +sT , + , 

-yo = ( v 2 ) z 2 + Y 2 z 0 + 


(48b) 


Any open branch arcs incident with Pj or P 2 are oriented towards these 
vertices, and their field relations are assumed to be included in eqs. 

(46a) and (47a), respectively, with v^ = v 2 = 0. In addition, the 
boundary conditions at Px and P 2 are [cf. eq. (2)]: ! 

_L _ | 

B i(I yi - l yi) + D izi = I (49a) 

B 2 (y£ ~ I y2> + D 2 z 2 = L 2 1 (49b) 

j 

*Tacitly included in this definition is the special case of a closed arc 
[e.g., in fig. 6(b) the circuit C if it consists of a single arc], in which 
case the "other nonpreceding" arc is the initial closed branch arc itself. 
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The solution of eqs. (46) -(49) may be conveniently divided into four 
steps. First, eqs. (46a) and (47b) are substituted into eq. (49a) to 
eliminate all yj and yj". This gives the desired relation 


2 i = v-2 z 2 + yo z o + y 

(50) 

y 2 = A - " 1 Bi £(v 2 ) T 

(51a) 

y 0 = A-^iV! 

(51b) 

y = a“ 1 [Li + Bi(J wi + £ x 2 )] 

(51c) 

A = Di - Bi (£ Ui + £ Y2) 

(51d) 


As a check of this result, note that if A is an open branch arc (vi - 0) 
and all are absent (v 2 = X 2 = Y 2 “ 0) > eq. (50) reduces to the formula 
for the terminal value of z for a tree [eq. (12) of ref. 1], It is noted 
also that eq. (50) is valid even if Pi is a singular vertex, since the 
inverse of Bi does not appear in eqs. (51). Upon reaching P 2 during the 
backward integration, both zq (the value of z at the I.P. of arc E) and 
Z 2 are known. The value zi (at the I.P. of arcs of the chains C^) is 
then calculated from eq. (50), whereupon the backward integration of eq. 
(12) on the chains can be made. 

Second, eq. (50) is substituted into eq. (47a) to eliminate z \ , and 
the result substituted into eq. (49b) to eliminate all y 2 . The result is 


B 2 y2 + D 2 Z 2 “ ^ 2 VZ 0 ~ 1-2 (52) 

where ^ 2 = D 2 - B 2 [£ u 2 4- (£ v 2 )y 2 ] (53a) 

u = (I v 2 )y 0 (53b) 

L 2 = L 2 + B 2 [l w 2 + (I v 2 )y] (53c) 


Note that if A is an open branch arc (u = 0) and E is absent (y* = 0) , 
eq. (52) gives z 2 » the terminal value of z, in a form similar to the 
formula for a tree. 

Third, the initial values of the field functions u, v, and w on the 
continuation arc E are obtained by substituting eq. (48a) into eq. (52) to 
eliminate y* and demanding that the resulting equation be an identity in 
zq and z 2 . The results are 

u 2 = -B 2 1 D 2 

4“ ^ 1 ^ 

w 2 = B 2 l ^- , 2 

+ 

v 2 - 
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(54a) 

(54b) 

(54c) 



Fourth, to obtain the initial values of y and A on arc E, eq. (50) is 
substituted into eq. (46b) to eliminate zj, and the result is substituted 
into eq. (48b) to eliminate yo. This result, also being an identity in 
zq and Z 2 , yields 

+ - - T 

Y2 = Yl + (vi) Po (55a) 

+ (v7) T m (55b) 

Note that if arc A is an open branch arc, then eqs. (54) and (55) imply 
that arc E is also an open branch arc, as it should be. For vf = yI = 

AT = 0 imply, through eqs. (51b), (53b), (54c), and (55), that v$ = y$ = 
A'J = 0. Also, it may be easily shown that eqs. (54a) and (55a) maintain 
the required symmetry of u and y. 

If P 2 is a singular vertex, eqs. (54) and (55) are replaced by 
formulas for the initial values of modified field functions. These may 
be obtained directly from eqs. (25), (34), (54), and (55) in a fashion 
similar to the derivation of eqs. (40). The results are 


u 2 - -(B 2 ot + D 2 3) _1 (D 2 ot - B 2 3) 

(56a) 

wt = (B 2 a + D23) _1 L 2 

(56b) 

vt = (B 2 ct + D23) -1 B 2 u 

(56c) 

-+ + T-+ 

y 2 = y 2 + u 3 v 2 

(56d) 

— Hh + T -f“ 

A 2 = A 2 + u gw 2 

(56e) 


where y 2 and A 2 are given by eqs. (55). As discussed on page 15, if E is 
a closed branch arc, or if E is an open branch arc with complete rigid 
body constraint prescribed at P 2 > then eqs. (56) can be simplified by 
choosing (3 = I [and hence by eq. (21) a = 0] . 

Now consider the case where the length h of one of the chains 
is very small. Assuming this chain to consist of a single arc, say F, 
and Pi to be a force-free vertex (i.e., Bj = I, Dj = Lj = 0), we expand 
the field functions on F in Laurent series about P]_. This gives the 
results 


u 2 = ~h 1 c 1 + 0(1) 

(57a) 

v 2 = h -1 c -1 + 0(1) 

(57b) 

p 2 = I + 0(h) 

(57c) 


where the matrix c may be evaluated at any point on r, and 0(1) means 
terms remaining finite as h 0. For small h, the term h -1 ^ 1 dominates 
0(1). Since the terms FT-’-c - * 1 cancel each other in the calculation of D 2 
[cf. eq. (53a)], it is clear that this calculation is subject to ever- 
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increasing round-off errors as h approaches zero. For this reason, it 
is necessary to treat the loop configuration [fig. 6(b)] as a special case. 

Loops: Initial values of the field functions on the continuation arc 
E of figure 6(b) can be obtained by applying the limiting process described 
above to eqs. (51)~(56). However, it is more direct to consider the loop 
configuration as a separate case and follow a procedure similar to that 
used to derive eqs. (52) -(56). Omitting the details, the results are the 
following. In place of eqs. (52) and (53) 


b 4 + 

Dz p - Bv A z 0 = L 


(58) 

where D = D 

~ B[ X u + l Y c + l v c + 

V T ] 

(59a) 

L = L 

+ B(£ w + l X“) 


(59b) 

In place of eqs. (54) and (55) 



+ 

u^ = - 
E 

-B -1 D 


(60a) 

+ 

W E = 

B _1 L 


(60b) 

+ 

V E = 

V I 


(60c) 

+ 

y e 

4 


(60d) 

4 = 

X I 


(60e) 

If P is a singular vertex, one obtains in place 

of eqs. (56) 


— + 
U E 

-(Ba + D3) -1 (Da - BB) 


(61a) 

-+ 

W E = 

(Ba + D3) _1 L 


(61b) 

_+ 

V E = 

(Ba + ]3 3 ) ~ 1 Bv~ 


(61c) 

— + 
y e 

y A + (v A )Tf3 4 


(61d) 

4 = 

a a + (v a )T65 e 


(61e) 


CALCULATION PROCEDURE 


As outlined previously (fig. 3), the calculation procedure essentially 
consists of three basic steps: 1) forward integration for the field 
functions u, w, v, y, and X,* 2) calculation of z at the final vertex, and 


^Unless otherwise specified, in this section u, w, v, y, and X will stand 
for u, w, v, y, and X, respectively, on singular arcs. 
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3) backward integration to obtain the response functions y and z. 


Forward Integration 

On regular (ordinary and continuation) arcs, the differential, 
equations (14) and (16) are integrated starting with the initial values 
given by eqs. (18) and (19) for ordinary arcs, (54) and (55) for 
continuation arcs exiting from a split circuit, and (60) for continuation 
arcs exiting from a loop. On singular arcs, the differential equations 
(14) and (16) written for modified (barred) variables are integrated 
starting with the initial values given by eqs. (40) for ordinary arcs, 

(45) for initial closed branch arcs, (56) for continuation arcs exiting 
from a split circuit, and (61) for continuation arcs exiting from a loop. 

As they are calculated, the primary field functions u, w, and v (u 
and w on open branch arcs) , which are used in the subsequent backward 
integration, are stored. At the terminal vertex of singular arcs, it is 
convenient to replace the modified functions u, w, and v by u, w, and v, 
respectively, using the inverse of eqs. (25). This is done so that the 
initial values of the field functions for all arcs, as well as the value 
of z at the final vertex (see below), may be computed the same way, 
regardless of whether preceding arcs are regular or singular. 

Either a Runge-Kutta or a predictor-corrector integration scheme may 
be used to integrate eqs. (14). Typically, the field functions have a 
narrow zone of rapid variation (i.e., a boundary layer) immediately 
following the initial vertex of each arc, but are otherwise slowly varying. 
In order to take advantage of this behavior, it is desirable that the 
integration scheme use a variable step size which is automatically 
controlled to meet a prescribed truncation error tolerance. Then the 
most efficient storage locations for u, w, and v (to obtain their best 
description for the least amount of storage) are the end-points of each 
integration step. In this way, the spacing of the storage points will 
vary with the rate of variation of the stored functions, the more rapid 
their variation, the closer together their storage points. In order to 
obtain the greatest possible accuracy, it is desirable to store, in 
addition to the primary field functions, their s-derivatives , which are 
avialable at no extra effort. Then a more accurate "double point" 
interpolation polynomial can be used to calculate intermediate values of 
these functions. The author has had good success with a fourth-order 
Runge-Kutta routine which proceeds in double steps, each of which is 
composed of two equally spaced basic Runge-Kutta steps.* In order to 
maintain fourth-order accuracy throughout, fourth-order interpolation of 
the primary field functions is made during the backward integration. For 
this purpose, these functions are stored at the end of each basic Runge- 
Kutta step, whereas their derivatives are stored only at the end of each 


*Quadrature of eqs. (16) for y and A is made over each double step by 
Simpson’s rule. 
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double step. Then Interpolation at points lying in the range of a double 
step (si < s < S 3 ) is made according to the formula 

f(s) = e 2 (e - 2 ) 2 f 2 - (0 - i)(e - 2 ) 2 [ (i + 2 e)f 1 + ehfjMM 

+ (0 - 1 ) 0 2 [ (5 - 20)f 3 + (0 - 2 )hf 3 ’ ]/4 + 0(h 5 ) (62) 

where 0 = (s - si)/h (63a) 

h = s 2 - si = S 3 - s 2 (63b) 

Here, fj_, f 2 , and f 3 are function values at the step points S]_, s 2 , and 
s 3 , respectively, and f^’ and f 3 f are s-derivatives of f at s\ and S3. 

With such a storage method, for any given problem the number of 
storage locations for the primary field functions is not known a priori. 

In order to avoid overflow of any fixed size storage area in core, and 
also to reduce the core storage requirement, a method of "dynamic" 
storage is used. In this method, only a relatively small buffer area 
for field function storage is set aside in core. When, during forward 
integration, this buffer fills up, it contents are automatically 
transferred to an external mass storage file. Because in a sequence of 
problems for which only load changes occur, u and v need not be recalcu- 
lated, it is convenient to actually use two external files, one for w 
and the other for u and v. Hence at the end of forward integration, in 
general, two external files consisting of a sequence of logical records, 
one for each buffer load, will have been created. During backward 
integration, these logical records are returned to the core buffer 
area (in reverse order) when the integration progresses out of the region 
of applicability of the current buffer contents. 


Calculation of z at Final Vertex 

There are three types of final vertices at which the calculation of 
z differs. The final vertex (i. e. , the terminal vertex of the final arc) 
may be: 1) the closure point of a split circuit, such as point P 2 in 
figure 6(a) with arc E absent, 2) the closure point of a loop, such as 
point P in figure 6(b) with arc E absent, or 3) incident with no closed 
branch arcs. 

Closure point of a split circuit .- In this case, arc A of figure 6(a) 
is necessarily an open branch arc, and z is given by eq. (52) with terms 
in y* and u omitted, viz. 

z 2 = D 2 1 L 2 (64) 

where D 2 and L 2 are given by eqs. (53a, c). 

Closure point of a loop .- In this case, arc A of figure 6(b) is also 

an open branch arc, and z is given by eq. (58) with terms in yt and vT 

Ei A 
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omitted, viz. 


Zp = D -1 L (64) 

A A 

where D and L are given by eqs. (59). 

No incident closed branch arc .- In this case, the formula for z may 
be obtained from eq. (50) considering all chains absent, or just as 
easily from eq. (58) considering all circuits C and arc E absent. The 
result is 

z - D” 1 ! = (D - B l u) -1 (L + B J w) (65) 

Note that this formula is more general than the corresponding formula, 
z = (D - Bu) -1 (L + Bw) , given in reference 1 for a tree. This is the 
case since the ordering rules for open branch arcs given in reference 1 
have been relaxed to the extent that the final vertex no longer need be 
a branch extremity. 


Backward Integration 

With the value of z at the final vertex as an initial condition, a 
backward integration of eq. (12) is performed over the graph, using z- 
continuity at interior vertices. Before integrating each closed branch 
arc, the stored values of w on that arc are replaced by w + vzq, where z q 
is the value of z at its I.P. Except for this replacement (and, of 
course, a similar one for w 1 ) there is no difference between the backward 
integration on open and closed branch arcs. Upon reaching point of a 
split circuit [fig. 6(a)], the value of zo for arcs of the chains is 
calculated from eq. (50). [In this case zq is z\ and should not be 
confused with the zq in eq. (50) , which is the known value of z at the 
I.P. of arc E.] Upon reaching point P of a loop [fig. 6(b)], the value 
of z 0 f° r arcs °f the circuits C is known since in this case, by z- 
continuity, z 0 = z p . From the values of z obtained by integration on 
regular arcs, y is calculated using eq. (7). 

On singular arcs, the backward integration is performed in terms of 
modified variables, i.e., eq. (12) written for barred variables is 
integrated. For this purpose, the value of z at the terminal vertex of 
a singular arc is replaced by z by first computing y there from eq. (7) 
(recall that u and w have been stored there) and then z from eq. (20b). 
After so doing, u and w at the terminal vertex are changed back to u and 
w through use of eqs. (25a, b), as required for the integration of the 
modified form of eq. (12).* From the values of z obtained by integration, 
y is computed from eq. (7) written for barred variables, and y and z are 


*Note that, because of the similarity of eqs. (25b) and (25c), the trans- 
formation from w + vzq to w + vzq on singular closed branch arcs is the 
same as that from w to w on singular open branch arcs. 
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then given by eqs. (23). 


SHELLS OF REVOLUTION 


A computer program employing the field method to obtain the linear 
elastic response of multicircuit ring-stiffened shells of revolution 
subject to general harmonic mechanical and thermal loads has been written. 
For this class of problems, the graph over which the boundary-value 
problem is defined represents the reference meridian of the shell. The 
differential equations (1) are eighth-order so that the response matrices 
y and z are 4-element column vectors. The equations are ordered so that 
y and z are the force and displacement vectors (fig. 7) 

y T - r(P, Q, S, Mj) (66a) 

z T = (S, n, v, x) (66b) 

where P, Q, and S are forces per unit of circumferential length referred 
to fixed axial, radial, and circumferential coordinate directions x, y, 
and <j> ; £, n, and v are the corresponding displacement components; M} is 
the meridional bending moment per unit of circumferential length; and x 
is the corresponding rotation. 

The matrices a, b, c, f, and g [eqs. (1)], as well as < = B _1 D and 
B~ 1 L [eq. (2)] for ring boundaries, are given in the appendix. Since 
this problem is self-adjoint, b, c, and k are symmetric and d = -a-*- [cf. 
eqs. (4)]. Corresponding shell and ring equations have been given 
previously (although not in precisely the same form) in reference 5. 


Example Solutions 

Solutions of the shell equations by the field method for three 
example multicircuit meridians are presented in this section. In each 
case, in order to have a quantitative check on the accuracy of the 
numerical solution, the structure was devised to have a plane of symmetry 
normal to the axis of revolution. (In practice, only one-half of such 
symmetrical structures would be modeled.) Thus, the accuracy of the 
solution can be checked simply by comparing the calculated response at 
symmetry points. 

Diagrams of the meridians of the three sample shell configurations, 
showing entry (output station) numbers,* branch point numbers (encircled), 
and harmonic pressure load amplitudes (underlined), are given in figure 8. 
In this figure, the arrow on each arc indicates the arc orientation, i.e., 
the direction in which the forward integration proceeds. Arc ordering is 


*Equally spaced output stations corresponding to intervening integral 
entry numbers are not shown. 
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indicated by the entry numbers in that during the forward integration the 
entry points are encountered in the order of increasing entry number. 

Note that the first and second examples are the same problem set up with 
different arc ordering and orientation. For this problem the circumferen- 
tial harmonic number of the pressure loading and response was n = 1000, 
whereas in the third example n = 10.* For both configurations the 
following homogeneous, isotropic wall properties were assumed: E = 10 7 , 
v = 0.25, t = 0.1. These properties and the dimensions and pressure 
amplitudes shown in figure 8 may be taken to be in any consistent set of 
units . 

Each of the problems has been run in single precision on the UNIVAC 
1108 computer (a 36-bit word machine) with no degradation of accuracy 
relative to the CDC 6600 computer (a 60-bit word machine). Therefore, it 
is concluded that the small numerical discrepancies which occur at 
symmetry points (see table I, which is reproduced from computer print-outs 
with an added column labeled "SYM. POINT” giving the entry number of the 
symmetry point for each output station) are due essentially to truncation 
error, which is controlled by an input truncation error tolerance, and 
that round-off errors are negligible. 

For each calculation the truncation error tolerance was set to 10~ , 
which means that at the end of each Runge-Kutta double step (see p. 21), 
in both forward and backward integrations, the relative difference between 
the values of each (scalar) dependent variable as calculated by the 
Runge-Kutta and Simpson rules (over the double step) was less than IQ -3 . 
Roughly speaking, this means that for each step these two calculations 
agreed to three significant digits. Since, as noted, round-off errors 
were negligible, one expects this degree of accuracy in the response 
functions (except in the immediate vicinity of their zeroes). Inspection 
of the calculated values of the components of the force and displacement 
vectors for each of the sample problems (table I) shows that this is 
indeed the case. 

Execution times (central processor times in seconds) and the number 
of double steps made during the forward and backward integrations for 
each problem are shown in the following table. 


Example Execution time Forward steps Backward steps 


a 15.2 (U-1108) 124 110 
b 14.8 (U-1108) 208 104 
c 12.6 (CDC-6600) 300 218 


*In concocting such hypothetical problems, it should be kept in mind that 
for a given shell the value of n cannot, within the scope of shell theory, 
be chosen arbitrarily large since nt/r must remain small compared to unity. 
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Comparison with Other Methods 


Conventional methods of numerical analysis of the response of shells 
of revolution suffer from several practical difficulties. In finite- 
difference methods one must choose in advance the finite-difference mesh 
points. The difficulty here is that too few mesh points may result in 
excessive truncation error, whereas too many mesh points may result in 
excessive round-off error. This difficulty is illustrated in figure 9, 
which has been taken from reference 6, giving the results of a bifurcation 
buckling analysis of an axially loaded cylindrical shell as a function 
of the number of finite-difference intervals on the shell half-length. 
Numerical convergence appears to be obtained with 27 intervals since the 
critical loads calculated with 27 and 28 intervals agree to four signif- 
icant digits. However, when the number of intervals is increased further, 
the growth of round-off errors causes divergence. 

In the superposition methods, which employ forward integration 
techniques, the truncation error is controlled automatically by the use 
of self-checking variable step-size integration routines. However, since 
the auxiliary solutions to be superimposed exhibit exponential growth 
characteristics, their superposition involves large round-off errors at 
points remote from the initial point of integration. To circumvent this 
problem, the meridian of the shell is artificially subdivided into a 
number of "suitably small" segments. One difficulty with this approach 
is that the subdivision requirements are known only very roughly, and 
depend on factors in addition to the properties of the structure itself, 
e.g., the circumferential wave number and, for eigenvalue problems, the 
eigenvalue shift. Also, some problems may not be treatable simply 
because more segments are required than allowed by computer storage 
limitations . 

The field method not only retains the advantage of automatic control 
of truncation error, but also produces stable initial value problems whose 
solutions are not sensitive to small numerical error introduced at each 
integration step. Consequently, no artificial subdivision of the shell 
meridian is required, leaving the user free from concerns of numerical 
convergence (truncation error) or ill-conditioning (round-off error). 

This is illustrated in figure 10, showing computation times for the cal- 
culation of the linear elastic response of a pressurized clamped cylinder 
by the field method and also by an efficient superposition method, known 
as the Zarghamee- Robins on method (ref. 7).* The bar graph [fig. 10(a)] 
compares the computational efficiencies of the two methods in terms of 
CDC 6600 CPU seconds for the axisymmetric response for several r/t ratios. 
In each case the field method is significantly more efficient. More 
important is the fact that, even with 33 segments, the superposition 
method failed to give meaningful results at r/t = 10^, whereas the field 
method experienced no difficulty without any subdivision.. Similar results 
are obtained when the circumferential wave number is varied, as shown in 
figure 10(b). 


*The author is indebted to Dr. Wendell B. Stephens of the NASA Langley 
Research Center for providing the results shown in figure 10. 
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To provide additional verification of the accuracy of the field 
method for closed branch shells, a comparison of the field and Zarghamee 
solutions for the shell depicted in figure 11 is presented.* (For 
comparisons of field and Zarghamee solutions for open branch shells, 
see reference 1.) The meridian of this shell consists of a single circuit 
with attached open branches. Only entry numbers at physical boundaries 
are shown in the figure. The edges at entries 15 and 53 are clamped, 
whereas the remaining four edges are force-free. The interior boundaries 
at entries 30, 68, and 76 have attached rings, whereas the remaining 
three interior boundaries are force-free. Each of the open branch 
segments (between boundaries) has a unit slant length, whereas each of 
the closed branch segments has a slant length of 2 units. The segments 
25-29 and 48-52 are cylindrical, whereas all others are conical at 45° 
from the axis. The following homogeneous, isotropic wall properties 
were assumed: E = 10 7 , v = 0.25, t = 0.1. The three rings are identical, 

with the following section properties: EA = 2. x 10 6 , EI X = 2.1667 x 10 4 , 

Ely = 1.6667 x 10% EI X y = 0, GJ = 5.6240 x 10 2 . For each ring, it is 
assumed that the section centroid lies on the shell reference (i.e., 
middle) surface. The only load consists of a uniform second harmonic 
(n = 2) unit pressure loading as shown in the figure. In table II 
force and displacement components calculated by each of the methods are 
compared at the boundary points of one-half of the symmetrical structure. 
For the results of the Zarghamee-Robinson method only those digits which 
differ from corresponding digits of the field method results are shown. t 
As the table shows, the agreement is very good, in spite of the use of a 
truncation error tolerance of only 10 -2 in each calculation. 


CONCLUDING REMARKS 


The field method of solution of general even-order linear boundary- 
value problems defined on one-dimensional domains containing circuits has 
been formulated. The method has been implemented in a computer program 
for the static elastic response of ring-stiffened orthotropic shells of 
revolution with multicircuit meridians. Example solutions have been 
presented which illustrate the accuracy and efficiency of the method. 

The field method has several advantages over more conventional 
numerical methods such as other numerical integration (superposition) 
and finite-difference methods. These are: 

(1) The field method is essentially free of numerical ill- 
conditioning problems. All calculations are made in single precision. 


*Since the Zarghamee-Robinson method has been extended to closed branch 
shells with only a single circuit (ref. 5), it could not be used to 
provide solution comparisons for the shells of figure 8. 

fin the table, an exact zero is indicated by 0., whereas 0.000 indicates 
a number which is zero to three decimal places. 
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even on 36-bit word computers such as the UNIVAC 1108, and there is no 
need for artificial subdivision of the domain of Integration. 

(2) Solution of large systems of algebraic equations is avoided. 

The method is developed in terms of matrices whose order equals the 
half -order of the boundary-value problem. 

(3) A priori choice of output or storage points for response 
variables is unnecessary since a uniformly dense set of storage points 
can be chosen automatically during execution. This is a particularly 
important feature when combined with iterative methods for nonlinear or 
eigenvalue problems, for which the response of each iteration is the 
input to the succeeding iteration. Consequently, there are no convergence 
problems related to an inadequate choice of response (or finite-difference) 
points. 

Because of these features, the field method is potentially a truly 
automatic method in that one can expect, with confidence, from computer 
programs based on it meaningful answers for every meaningful problem 
for which the program is designed. 
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APPENDIX 


MATRICES FOR AXISYMMETRIC SHELLS AND RINGS 


A 

E 


e x ,e y 


F x’ F y’ F < 


GJ 

I X » ly » 
k 

Li,L 2 




N x ,N y ,^ 


Xi,X 2 ,X 3 


o <°> n (0) 

L>1 ,0 2 


n (1 > n (1) 

©1 ,© 2 


X ij 




Symbols 


ring section area 

ring elastic modulus 

ring centroidal eccentricities 

harmonic amplitudes of ring force loads per unit of 
circumferential length 

ring torsional stiffness 

ring section moments of inertia 

ring stiffness matrix 

harmonic amplitudes of shell moment loads per unit of 
surface area in meridional and circumferential directions 

ring load vectors 

harmonic amplitudes of ring moment loads per unit of 
circumferential length 

harmonic amplitudes of shell force loads per unit of 
surface area in meridional, circumferential, and 
normal directions 

ring eccentricity matrix 

harmonic amplitudes of meridional and circumferential 
thermal force loads 

harmonic amplitudes of meridional and circumferential 
thermal moment loads 

harmonic amplitude of ring free thermal strain 
orthotropic shell wall normal stiffness coefficients 
orthotropic shell wall shear stiffness coefficients 
ring centroidal radius 
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Shell 


For elastic orthotropic shells of revolution subjected to symmetric 
n-th harmonic loads, the coefficient and load matrices of eqs. (1) are 
given below. For the definitions of the stiffness (A^j, ^ij) and load 
[X^, L^, variables in these matrices, see reference 5. 
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ft IH ft f H-* H |}-i 


n2 A 2 3 r ’/R2 

-(A 13 +n 2 A 23 /R 2 ) (r/R 2 ) 
-n(A 13 +A 23 /R 2 )(r/R 2 ) 
(1-A 23 /R 2 )rr T 


n2 A 23 r 1 2 /r 

n(r/R 2 -ui 2 /r) 

n 2 A 2 4r'/r 

-(A 1 3 +n 2 A 23 /R 2 )r' 

nr’ 

~ (Am+n 2 A 2 4/R 2 ) 

-n(Ai 3 +A 23 /R 2 )r f 

r’ 

~n(Ai 4+A 2 4/R 2 ) 

-r 2 /R 2 -A 23 r ? 2 

ny 12 

" X 24 r ' 


n 2 (n 2 A 22 r ’ 2 +yn)/r 2 -n 2 (A 12 +n 2 A 22 /R 2 )r '/r 

A i i+2n 2 A ^ 2 /R 2 +n l+ A 22 /R 2 2 

Symmetric 


-n 3 (A 12 +A 22 /R 2 )r'/r 

n[A 11 +(n 2 +l)A 12 /R 2 'hi 2 A 22 /R 2 2 ] 

n 2 (A n +2A 12 /R 2 +A 22 /R 2 2 ) 


“n 2 (A 22 r l2 +u u )/r 

(A 12 +n 2 A 22 /R 2 )r’ 

n(Ai 2 +A 22 /R 2 )r 1 

A 22 r ,2 +n 2 un 


f = 


A 33 (r/R 2 ) 2 A 33 r T (r/R 2 ) 

*33r 12 

L Symmetric 


0 

0 

P22 


A 3 4 (r/R 2 ) 

^34 r ? 

0 

X 44 J 


f rt-(r/R 2 )X 1 -h: , (X3-nL 1 /r)]+n 2 (r , /r)[02 (1) -X230i (O) -X2401 (1> I 

-r[r , X 1 +(r/R 2 )(X3-nL 1 /r)]+[X 1 3+n 2 X23/R2]0i (O) +[Xm+n 2 X24/R 2 ]0i (1) -e2 (O) -n 2 02 (1) /R2 

-rX2+(r/R2)L 1 +n[(X 1 3+X 2 3/R 2 )0 1 (O) +(X 1I( +X 24 /R 2 )0 1 (1) -0 2 (O) -0 2 (l:) /R2] 


lO 




-rL 2 +r ' I X 2 3©i ^°^+X 2 40 i 
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, the coefficient and load 
(ref. 5) 


n 3 EI xy /p 2 -n 2 (Ely+GJ) / p 

n(EA-Hi 2 EI x /p 2 ) -n 2 EI xy /p 
n 2 ( EA+EI X / p 2 ) -nEI xy / p 

EI y +n 2 GJ 


“ N x> P N ^ 
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TABLE I.- CALCULATED FORCE AND DISPLACEMENT AMPLITUDES 
(a) Split Circuit in Loop 


SYM. 


POINT 

FwT pV 

S 

0 

0 

CAP F 

M] 

XI 

e t a 

V 

CHI 

31 

1 

ft . ftono 

7.1078+02 

-1 ,3701-03 

2.2391.02 

-1.51 40+00 

9,6569-00 

-1.0752-03 

-3.5032-06 

-1.0213-07 

30 

•? 

l . 1 ) 07 + 00 

p.ltso+o? 

9, ?F?P +0 1 

-1 .0557 + 02 

2 . 1 9a« + 00 

4.5550-04 

-1.7270-03 

-4.3766-05 

-5.2469-04 

29 

7 

9.7?| 4 + ft ft 

7.5997*0 1 

6 # 5775+0i 

-1 ,8035 + Op 

.6.9776*00 

1 .2166-04 

-1,0245-03 

4,0172-05 

3.8299-04 

28 

4 

9.2214 +nn 

9.7 33F + 0 1 

2 • 8 1 06 + 0 i 

-4 .8407 + 0) 

6,6247+00 

1.2168-04 

-1.0245-03 

4.0172-05 

3.6299-04 

27 

f 

7,O06A*.<)O 

A ,8749 + 01 

4,7739+01 

-5.0880+01 

-1,71 42+00 

-1,6699-04 

-1.7255-04 

4, J005-06 

6,3368.04 

26 

ft 

7 . 7 97 p * ft ft 

• 4 ,*165-03 

4.27U+01 

-) .4928-6, 

-4,0186-0! 

*3.4?59*04 

-6,65) 0-06 

9,6196-05 

4.6671-07 

25 

7 

4 .F77ft+on 

•4.6707+01 

4.7728+01 

5,0527+0 1 

-1.7117+00 

-1.6916-04 

1 .7260-04 

4,2157-06 

-8.3360-04 

24 

a 

F.76.Jft'‘- n ft 

-9,7105+0) 

2 . 8028+01 

4.8049+0) 

8 , 61 13+00 

1.2169-04 

1.0242-03 

3,6445-05 

-3.6341-04 

23 

0 

F . 7ft3o+oo 

-6, 1335+01 

5.7634+nl 

-1.3196+02 

-) .3600 + 01 

1,2168-04 

-1.0245-03 

4,0172-05 

3.8299-04 

22 

in 

0 . 1 484+ftft 

1 .2580 + 0) 

-1.3814+01 

-1 .0727 + 02 

2.82)9+00 

5,1937-04 

-7,3955-05 

6.2660*05 

-8.5601-04 

21 

1 t 

ft.Q33A+00 

A .4554-03 

-6.7498+01 

-1 ,8100-0! 

-2,3092-01 

6,2279-04 

-1 .6187-07 

1.4441-04 

-1.2323-06 

20 

1 9 

7 , 7 j 97 + 00 

-1 .2643+0 1 

-1.7855+01 

) ,0697 + 02 

2.6245+00 

5.2031-04 

7,3276-05 

8.1973*05 

8.5764-04 

19 

t 7 

A ,F040+ftft 

0 . ) -370+01 

5.7642+0! 

1 ,3l«3 + 02 

-1.3607+01 

1,2169-04 

1.0242*03 

3,6445*05 

-3.8341-04 

18 

1 4 

A . SO 4A + 00 

-3.5834+Oj 

8.5670+01 

) .7986 + 0? 

-6.9964+00 

1,21 69-04 

1,0242-03 

3.6445*05 

-3.8341-04 

17 

IF 

O.(Sl*i7 + 00 

-7.3687+02 

9,2610+0 i 

1 .0577*0? 

2,2071+00 

4,5716-04 

1,7289-03 

-4.6603*05 

5.2704*04 

16 

1 * 

1 . 0 / 20 + 0 ) 

-7.1151+02 

3.5245-04 

6.3405-03 

-1 .5224+00 

-1.6873-08 

) ,07?9*03 

-5,5312*05 

9,4098-08 

15 

t^ 

1 . 1 «37 + 0 t 

-7.3657+02 

-‘ 1. 2606 + 0 ) 

-) .0576*0? 

2,2071+00 

-4.5709-04 

1.7267-03 

-4,6602*05 

•5.2772*04 

14 

1 * 

1 .994« + 0 t 

-7,581 5+0) 

-8 , 5660+0 l 

-1.7967+0? 

-6,9971+00 

*1.2194*04 

1.0242-03 

3,6435*08 

3.8331-04 

13 

i 0 

1 . 794A+0 i 

ft . 1 280 + 0 J 

-5.7645+0) 

-1.3165+0? 

- l , 3606*0 l 

-1,2194-04 

1.0242-03’ 

3,6435*05 

3.6331*94 

12 

20 

1 .3737+01 

- ) ,2640+0 1 

1.3653+0) 

-1 .0700 + 0? 

2.8247+00 

-5.1967-04 

7,3359-05 

6,2025*05 

*8.5639*04 

11 

21 

1 , 4*i 1 « + 0 1 

7.758 1-03 

6.7499+01 

, .6648-0) 

-2,3096-01 

-6.2281-04 

- 3 , 211 6 -O 7 

1.4442*04 

3.4690*07 

10 

2 ? 

1 • F 3 0 4 ♦ 0 i 

1 .2580 + 0] 

1.3814+0) 

1 .0725+0? 

2,6220+00 

•5, 196 1-04 

-7 , 39)0-05 

§,267?»0B 

6.5695*04 

9 

?7 

) , 0089+C | 

-6.1325+0) 

-5.7631+01 

1 .31 94 + 0? 

-1 .3601+01 

-1,2163-04 

-1.0244-03 

4,9180*05 

-3.8314*04 

8 

7 A 

1 . ftOftq+O) 

-0.7099+0) 

-2.8ft96+0i 

-4 . 6028+0 1 

6 , 6106+00 

-1.2194-04 

1.0242*03 

3,6435*05 

3.6331*04 

7 


1 .FA7F + 0 , 

-4.67 ) O+ft ] 

-4.7727+01 

-5.0532+0, 

-1 .7117*00 

1, 6868-04 

1,7287-04 

472rST*93 

6^3406*94 

6 

2 * 

1 • 7 ft60 + ft 1 

-7.6580-03 

-4.27) 1 +01 

9.8543-0? 

-4.0153-0) 

3.4255-04 

2.7275-07 

9.6193-05 

4,6755*07 

5 

27 

1 . 644F+0) 

4.6743+01 

-4.7734+01 

5.0670+0) 

-1,7)41+00 

1.6939-04 

- 1 • 7 231-04 

4,3036*05 

— 6 . 3 3 1 6*04 

4 

7* 

1 .Q?3l + 0) 

9.7338+01 

-2.8107+01 

4,6421 +0) 

6.6248*00 

-1,2163-04 

-1.0244*03" 

479189*05 

831 

3 

?0 

1 .9231 *0) 

7.601 3 + 0) 

-5.5737+01 

) .6036+0? 

-6,9766*00 

-1.2163-04 

-1.0244-03 

4,0169*05 

-3.6314*04 

2 

1 ft 

9.0.14 ) + 0 ) 

7.3659+02 

-9.2530+01 

) .056, +0? 

2,1951+00 

— 4 . 5563—04 

-1,7273-03 

•4,3629*05 

5,2474*04 

1 

3i 

7.1 457+0 1 

7 . 1 080+02 

-1.058 1 -03 

-4 , 8722-03 

-1,5135+00 

9.8569-09 

-1.0752-03 

-3.5032*06 

*I.0Z13i97 
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TABLE I.- CALCULATED FORCE AND DISPLACEMENT AMPLITUDES - Continued 


(b) Split Circuit in Split Circuit 


STM. 


POINT 

21 

ENT*Y 

1 

2 

S 

0.0000 

P 

9.7331+01 

0 

2.6l07+0i 

CAP 8 

-4.8403+01 

6 

Ml 

•6246+00 

XI 

1.2166 

*1,6900 

04 

ET* 

•1.0245-03 

V 

4.0175-05 

4+3O68»0f 

CHI 

3. 8306-04 

22 

7.8540*01 

4*6746+0i 

4.7740+0i 

•5.0872+0i 

•1 

.71*4+00 

04 

•1*7^2 4 9*04 

ft.336^04 

23 

3 

1,5708+00 

•2.6678-03 

4, 2715+01 

•1 ,5l67-0i 

•4 

, 0 1 69-0 j 

•3.4J47 

04 

-3,2087-01 

9.6169-05 

1,5466-07 

24 

4 

2,3562+00 

•4 , 6705401 

4,7723+01 

5.0520+01 

-1 

»7i 15 + 00 

•1.8911 

04 

1.7J82-04 

4.2146-05 

-6,3361-04 

25 

6 

3.1416+00 

•9.7095+01 

2.6024+01 

4,6029+01 

6 

•6106*00 

1.2193 

04 

1.0242-03 

3.844X-05 

-3,6344-64 

16 

ft 

3.1416+00 

•6.1340+01 

5 , 763ft+0 1 

-1.3196+02 

-1 

.3600+01 

1.2166 

04 

•1.9245-03 

4.0175*05 

3.8308-04 

17 

7 

3.0270*00 

1 .2578+01 

•1.38l3+0i 

•1.0728+02 

2 

.6221*00 

5.1931 

04 

•7.4052*05 

6,2660-05 

-6,5771-04 

18 

ft 

4,7124+00 

8.9391-03 

-6,7505+01 

• l ,9449-0 1 

•2 

.3|37-0i 

6.2274 

04 

•2, 3260.07 

1,4440-04 

•1 ,34*2-06 

19 

0 

3.4976+00 

•1.2639+01 

-1.3653+01 

1,0699+02 

2 

•6249+00 

5,2036 

04 

7.3204-05 

6.2018-05 

6,5755.04 

20 

10 

6.2832+00 

6.1275+01 

5.7645+01 

1.3184+02 

•1 

,3607+Ot 

1.2193 

04 

1*0242*03 

3+6441-05 

-3,6344-04 

15 

It 

6,2832+00 

-3.5979+0 l 

-8,57l9+0i 

1,8037*02 

6 

,9733+00 

1.*1H 

04 

•1,0245-03 

4.0i75-05 

3.6306-04 

14 

I? 

7.3939+00 

•2.3657+Og 

-9.25l9+0i 

1 ,0563*02 

•2 

•1964+00 

4,5534 

04 

•1 .7267*03 

*4+3605-05 

-6,2606-04 

13 

13 

6.504 6 +00 

•3,1080+02 

3,2243-03 

* .9967*03 

1 

.5 1 30+00 

1.5 6 05 

0? 

•1.074ft. 03 

-3.4381-Ofl 

-t, 2270-06 

12 

14 

9,6l53*00 

*2.3661+02 

9,2539+01 

-1,0549+02 

-2 

« 1949+00 

*4.5527 

04 

*1 +7267-03 

-4+3587-05 

5,2457-04 

11 

13 

l,0726+0i 

•3,5996+01 

8.5739*01 

*1.8034*02 

6 

.9761+00 

-1-2167 

04 

•1 .0245*03 

4.0154*05 

-3,6304-04 

6 

13 

1,0726*01 

6.1337+01 

5.7637+01 

-1.3195*02 

1 

•380i+0i 

•1.2167 

04 

•1 +0245-03 

4. 0154*05 

-3.6304-04 

7 

\y 

1 . l5i l+Oj 

•1 .2576+01 

*1,3812+01 

-1 ,0726+02 

•2 

•6221+00 

•5,1936 

04 

•7,3972*05 

6, 2671*05 

6,5803-04 

8 

13 

1.2297*01 

•8.491 1»0J 

-6,7503*01 

-1,9469-01 

2 

•3» 11-oi 

•6.2279 

04 

•2+0201-0? 

1.4440-04 

1.2462-06 

9 

to 

1,3062*01 

l,3666+0i 

1.2644+01 

,1,3855*0i 

1,0697+02 

• 2 

,6244+00 

•5,2029 

04 

7,3286-05 

6.1972-05 

-6,5764-04 

10 

20 

-6.1267+01 

5,7640+0i 

1.3183+02 

1 

,3606*01 

•1.2186 

04 

1,0242-03 

3.8447-05 

3,8336-04 ' 

1 

2t 

l,3 fl 6a+oi 

•9.7339+01 

2,6112+01 

-8.6*02+01 

•6 

,6254+00 

-I.2167 

04 

•1,0245-0? 

4.0154-05 

-3.69Q4-04 

2 

2? 

i , 4653+0i 
l,3436+0i 

•4.675?+0i 

4,7742+01 

■O.O666+O1 

1 

• 7 1 45 + OO 

1.6904 

04 

•1, 7250-04 

4.2*95*05 

-6,3377-04 

3 

23 

3,6420-03 

4, 2 7i9+0i 

-1 .3743-01 

4 

,02l0-0f 

3.4271 

04 

3. 4663*08 

9.6I?!fi05 

-3,4637-07 

4 

24 

1. $224*01 
1, 7009+Qi 
i, 7009+Qi 

4, 6705+01 

4,772?+0i 

5,0523+01 

1 

•7 1 16+OO 

1.6909 

04 

1*7262-04 

4.2161-05 

6,33*2-04 

5 

23 

9,7i07+0i 

2,6028+0i 

4.8045*01 

-6 

•6j 15+00 

•1.2186 

04 

1.0242*03 

3.6447-05 

3.6336-04 

30 

2ft 

3,5639+0i 

8,5669+0i 

1.7988+02 

6 

.9947+00 

•1.2118 

04 

1.02X2-03 

3*6447-05 

3,8S31i64 _ 

29 

27 

l.«120+0i 

2.3684+02 

9,2612+01 

1.0577+02 

•2 

•2070+00 

•4.S7J4 

04 

1. 7269-03 

•4.6599-05 

-5.2703-04 

28 

23 

1,9231+01 

3.1152+02 

2,4131-03 

1,5466-02 

1 

•5221+00 

4,6501 

06 

1.0730-03 

•5.5293*06 

6,8681-08 

27 

20 

2,0341+01 

2.3683+02 

•9,2606+01 

•1,0576+02 

•2 

,2070+00 

4,5712 

04 

i~7288i03 

»4;O60Si05 

H, 2763-04"' 

26 

30 

2,1452+01 

3.5620+01 

•8, 5666+01 

-1.7967+02 

6 

,9966+00 

1.2193 

04 

1.0242-03 

3.6441-05 

-3.6344-04 


•FIN 


CO 

at 



CO 

cn 


TABLE I.- CALCULATED FORCE AND DISPLACEMENT AMPLITUDES - Concluded 
(c) Loop in Split Circuit 

SYM. 


POINT 

ENTRY 

S 

P 

Q 

CAP S 

Ml 

XI 

ETA 

V 

CHI 

37 

1 

0.E+00 

4 , 69666*00 

9.23466*01 

•1 ,6702E*0l 

1.1 177E*00 

-7.7628E-07 

2,19586-05 

-2.69256-05 

-6. 0202E-06 

36 

2 

7.8540E-01 

-6.8064E*0l 

6.80636*01 

-1.26286*01 

-3,20556-0 1 

6.5584E-05 

1.28146-04 

-2.6667E-05 

T.5061E-05 

35 

3 

1.57086*00 

-9,94476*01 

-3.50846-01 

-1.36096*00 

-3.5185E-03 

5.8907E-07 

1,43486-04 

-2.5374E-05 

-5,04226-05 

34 

4 

2.35626*00 

-7.05666*01 

-7 . 0326E* 0 1 

6, 8517E*00 

1.15036-02 

-6.3112E-05 

9,63996-05 

-2.8753E-05 

4.9203E-07 

33,38 

5 

3. 14166*00 

6.3655E-02 

-9,9387E*01 

1.00666*01 

-5.79616-02 

-9.15146-05 

2.7534E-05 

-2.95886-05 

-1,36366-05 

39 

6 

3.9270E*00 

7,0783E*01 

-7 . 0620E*01 

6.7149E*00 

7 , 2900E-02 

-5.2810E-05 

-4.0263E-05 

-2,64686-03 

5.41306-06 

40 

7 

4.7124E*00 

9,98946*01 

-5. 7729E-0 1 

1.41616*00 

-4.38476-02 

8.4095E-06 

-1.0578E-04 

-1.5B36E-05 

7.5130E-05 

41 

a 

5,4978E*00 

6.80046*01 

6, 8353E* 0 1 

-1.4485E*00 

-3.4530E-01 

7.4569E-05 

-9.0601E-05 

-1.6604E-05 

-9.59606-05 

42 

9 

6.28326*00 

-4.24826*00 

9. 3188E*01 

-4, 1802E*00 

1,0463E*00 

-7,76286-07 

2.1958E-05 

-2.69256-05 

-6.02026-06 

43 

10 

6.2832E*00 

-8.94236*00 

8.38226-01 

1.2518E*0l 

-7.0639E-02 

-7.7628E-07 

2.1958E-05 

-2.6925E-05 

-6.0202E-06 

44 

li 

6.7832E*00 

-1 .22636*01 

-2.6804E-01 

1 .3231E+00 

5.4458E-02 

-3.50926-06 

2.34046-05 

-1.9590E-05 

1.15636-05 

45 

1? 

7 • 28326*00 

-1.02076*01 

-1.4214E*00 

-1.0168E*01 

-3.4806E-01 

-6.49176-06 

2.2499E-05 

-2.7311E-05 

-4.5062E-05 

27 

13 

7 ■ 28326*00 

-6.23036*00 

9,1580E*01 

-1.3722E*01 

-1,3790E*00 

-6.49176-06 

2.24996-05 

-2.731 IE-05 

-4.5062E-05 

26 

14 

8 • 0686E*00 

6.5782E*0l 

6.54376*01 

-9.72676*00 

6. 51126-01 

-1.O130E-04 

1.30B9E-04 

-2.8201E-05 

-5.3976E-06 

25 

15 

8.85406*00 

9 . 7632E*0 1 

-4,1 454E *00 

-7,6982E*00 

-1.0968E*00 

-5. 04716-09 

6.81426-05 

-4.2072E-05 

9.56626-08 

24 

16 

8.85*0E*00 

9.29196*01 

-4.14?5E*00 

-7,6931E*00 

1.1414E*00 

-5.04T1E-09 

6.8142E-05 

-4.2072E-05 

9.5662E-08 

23 

IT 

9,6394E*00 

6.85796*01 

6,7959E*01 

-1 ,43816*01 

-3.3161E-01 

1.23596-04 

-2.33276-05 

-3.1728E-05 

1.38626-04 

22 

18 

1.0425E*01 

-5.52566-01 

1.0015E*02 

-1.11506*01 

-1 .0458E-01 

1.56266-04 

3.1952E-05 

-2.29236-05 

-7.52976-05 

21 

19 

1 • 1210E*01 

-7.07416*01 

7 • 07?4E*01 

•4,47666*00 

6.4669E-02 

7.6222E-05 

8.2322E-05 

-2.9616E-05 

-4.35296-05 

20 

20 

1.19R6E*0l 

-9.9182E*0l 

-4.3515E-04 

1.71156-03 

1.0815E-02 

1.6584E-07 

9.9575E-05 

-3.3616E-03 

-9,66026-07 

19 

21 

1 »2781E*01 

-7,07406*01 

-7.0725E*01 

4.4706E*00 

6.4498E-02 

-7.5464E-05 

8.20116-05 

-2.9659E-05 

4.20366-05 

18 

?2 

1.3566E*0l 

-5.5293E-01 

-1.00166*02 

1,1140E*01 

-1.0440E-01 

-1.5548E-04 

3.2239E-05 

-2.2931E-05 

7,53286-05 

17 

23 

1.4352E*01 

6,85746*01 

-6.7957E*01 

1.4302E*O1 

-3.3235E-01 

-1.2337C-04 

-2,32526-05 

-3.1759E-05 

-1,38376-04 

16 

24 

1.51 376*0 1 

9.29126*01 

4, 1461 E*00 

7,6961E*00 

1.1415E*00 

-5.0471E-09 

6,01426-05 

-4.2072E-05 

9.56626-08 

15 

25 

1.5137E*01 

9.7625E*01 

4.1471E*00 

7.6814E*00 

-1.09756*00 

-5.0471E-09 

6.8142E-05 

-4.2072E-05 

9.5662E-08 

14 

?6 

1.5923E*01 

6.57846*01 

-6.5*38E*01 

9,7150E*00 

6.5092E-01 

1.0130E-04 

1.30876-04 

-2.8186E-05 

5.20166-06 

13 

27 

1.67086*01 

-6,2264E*00 

-9.15866*01 

1.3716E*01 

-1.3777E*00 

6.52046-06 

2.250IE-05 

-2.7310E-05 

4.5249C-05 

32 

28 

1.670BE*01 

-3.98636*00 

-9.2994E*01 

3.53776*00 

1,03036*00 

-6.4917E-06 

2.2499E-05 

-2.7311E-05 

-4.30626-05 

31 

?9 

1.74936*01 

6,86626*01 

-6.8926E*0l 

1.38986-01 

-2.17766-01 

-5.5693E-05 

-8.26S3E-05 

-1.7425E-05 

9.79276-05 

30 

30 

1.82796*01 

1.0085E*02 

6.5333E-04 

4.94116-03 

-7.83226-02 

-9.3645E-08 

-1.34826-04 

-7.6005E-06 

-2.99216-07 

29 

31 

1 .90646*01 

6,8663E*01 

6,8927E*01 

-1,58036-01 

-2.1764E-01 

5.5372E-05 

-8.2366E-05 

-1 , 7460E-05 

-9.82196-05 

28 

32 

1,98506*01 

-3 * 98066*00 

9.300RE.01 

-3.55396*00 

1.0290E*00 

6.5204E-06 

2.2501E-05 

-2.73106-03 

4.5249E-05 

5 

33 

l,9850E*0i 

6.1557E-02 

9.93816*01 

-1.00886*01 

-5.5735E-02 

9.2364E-05 

2.76256-05 

-2.9606E-05 

1.1B19E-05 

4 

34 

2, 0635E*0l 

-7.05686*01 

7,03296*01 

-6,84656*00 

1.3015E-02 

6.34356-05 

9.6605E-05 

-2.8743E-05 

-1.54616-06 

3 

35 

2.1420E*01 

-9,9445E*0l 

3,50046-01 

1,35946*00 

-3.4288E-03 

-5.5277E-07 

1.43426-04 

-2.5400E-05 

4.97936-03 

2 

36 

2»2206E*0 1 

-6.8067E*01 

-6,8067E*0l 

1.26126*01 

-3.2040E-01 

-6.53586-05 

1.2791E-04 

-2.6673E-03 

-7.54186-05 

1 

37 

2.29916*01 

4.69?3E*00 

-9,2356E*01 

1 , 6693E*01 

1,1 164E*00 

8.0537E-07 

2.19796-05 

-2.6903E-05 

5.92366-06 

5 

30 

2.2991E*01 

-5.37826-02 

-9,9373E*01 

1.0072E*01 

5.54666-02 

9.2364E-05 

2.7625E-05 

-2,96066*05 

1.18196-05 

6 

39 

2,3777E*01 

-7.07776*01 

-7.0616E*01 

6.6937E*00 

-7.2456E-02 

5.31556-05 

-4.07296-05 

-2.6422E-0S 

-5,99616-06 

7 

40 

2.4562E+01 

-9.9894E*01 

-5.7605E-01 

1,41746*00 

4, 3948E-02 

-8.3064E-06 

-1.0598E-04 

-1.5807E-05 

-7.48896-05 

8 

41 

2.5347E*01 

-6.8004E*01 

6.83546*01 

-1 .4532E*00 

3.4513E-01 

-7.4400E-05 

-9.0440E-05 

-1.66256-05 

9.6247E-05 

9 

42 

2.61 33E*0l 

4.2462E*00 

9.3193E*01 

-4,1786E*00 

-1.0460E*00 

8.0537E-07 

2,19796-05 

-2.6903E-05 

5.92366-06 

10 

43 

2.61336*01 

8. 93R6E*00 

8.3679E-01 

1,25156*01 

7 . 04826-02 

8.0537E-07 

2.19796-05 

-2.6903E-05 

5.92366-06 

11 

44 

2.66336*01 

1.2260E*0l 

-2.6873E-01 

1,3275E*00 

-5.4190E-02 

3,53846-06 

2.3380E-05 

-1.9581E-05 

-1 . 1613E-05 

12 

45 

2.7133E*01 

1.02076*01 

«1,4220E*00 

•1 , 0162E*01 

3, 4866E-0 1 

6.5204E-06 

2.2501E-05 

-2.73106*05 

4,52496-05 



TABLE II.- COMPARISON OF CALCULATED FORCE AND DISPLACEMENT AMPLITUDES 


Entry 3 

Method 5 

PxlO 

QxlO 

Sxl 0 

MixlO 2 

exio 6 

nxlQ 5 

vxlO 6 

X x l 0 5 

1 

field 

-2.738 

1.381 

2.751 

4.803 

0.009 

1.922 

-2.730 

0.001 


Z-R 

7 

79 


6 

0 



0 

9 

field 

-4.687 

3.628 

5.006 

-7.219 

-9.001 

1.164 

-5.893 

-2.888 


Z-R 



12 

8 

9 

5 

5 

7 

10 

field 

0. 

0. 

0. 

0. 

101.84 

9.889 

-5.586 

16.342 


Z-R 





6 

91 

8 

4 

14 

field 

-8.208 

-2.712 

7.996 

-31.827 

1.596 

0.023 

-3.560 

6.092 


Z-R 


1 

7 

5 



1 

3 

15 

field 

13.893 

9.796 

-12.263 

-10.318 

0.000 

0.000 

0.000 

0.000 


Z-R 


5 

8 

23 

0. 

0. 

0. 

0. 

19 

field 

12.800 

7.209 

-16.619 

23.377 

1.596 

0.023 

-3.560 

6.092 


Z-R 

799 


25 

1 



1 

3 

20 

field 

4.592 

4.498 

-8.622 

-8.450 

1.596 

0.023 

-3.560 

6.092 


Z-R 

1 

7 

8 

4 



1 

3 

24 

field 

3.411 

3.401 

-9.955 

-7.943 

-9.001 

1.164 

-5.893 

-2.888 


Z-R 

09 

399 

61 

2 

9 

5 

5 

7 

25 

field 

0. 

0. 

0. 

0. 

-9.932 

8.847 

-9.503 

-8.929 


Z-R 





9 


7 


29 

field 

1.345 

-5.913 

8.183 

20.731 

-9.001 

1.164 

-5.893 

-2.888 


Z-R 


4 

2 

2 

9 

5 

5 

7 

30 

field 

0.072 

1.082 

3.158 

5.427 

-9.001 

1.164 

-5.893 

-2.888 


Z-R 

0 

1 

9 

32 

9 

5 

5 

7 

38 

field 

0.888 

0.000 

-0.003 

1.701 

0.000 

0.239 

-4.861 

0.001 


Z-R 

7 


0 

2 



2 

0 


“Entry numbers refer to output stations shown in figure 11. 

L 

Tor the Zarghamee-Robinson method, only digits which differ from corresponding digits of the 
field method are shown. 





a 



Note: Point a may be the I.P. for arcs 1, 5, and 7. 

Point b (arc 2) is the I.P. for arcs 2 and 3. 

Point b (arc 4) is the I.P. for arc 4. 

Point c is the I.P. for arc 6. 


Figure 2.- Illustration of initial points and arc ordering and orientation 
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I 



Figure 3.- Flow chart of field method 
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/ 
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Figure 4.- Illustration of an ordinary arc 



Admitted 


Not admitted 



Figure 5.- Examples of graphs admitted by limitation of graph generality. 

The heavy data depict vertices joined by three distinct chains 
one vertex lies in a circuit not containing the other. 
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(b) Split circuit in split circuit 



(c) Loop in split circuit 
Figure 8.- Example multicircuit meridians 
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(a) Effect of radius-to-thickness ratio 


Figure 10.- Computation times of field and Zarg 
for pressurized cylindrical shells 


"4 



(b) Effect of circumferential wave number 


i-Robinsot] methods 






Figure II.- Closed branch shell used for comparison of 
field and Zarghamee-Robinson solutions 


NASA- Langley, 1975 CR-2535 




